// Zeta function
// Copyright 2020 Glenn McIntosh
// licensed under the GNU General Public Licence version 3
#pragma once
#include "gamma.h"
#include <cmath>
#include <complex>

// uses modified Joseph Ser series
// N is number of terms in series
template<int N>
class Zeta
{
public:
	Zeta()
	{
		// calculate prescaled binomial coefficient sums
		for (int n = 0; n <= N; ++n)
		{
			int64_t b = scale;
			for (int k = 0; k <= n; ++k)
			{
				binom[k] /= 2;
				binom[k] += b;
				b -= binom[k];
			}
		}
	}
	template<typename T> T operator()(T s)
	{
		bool invert = std::real(s) < 0.;
		if (invert)
			s = 1.-s;

		// calculation
		T result{};
		for (int k = N; k >= 0; --k)
			result += double(binom[k]) / exp(log(k+1)*s);
		result /= scale*2;

		// scale
		result /= (1.-exp(-s * log(2))*2.);

		if (invert)
		{
			constexpr double tau = ::acos(-1)*2;
			result *= gamma(s) * 2./exp(log(tau)*s)*cos(s*(tau/4));
		}

		return result;
	}
private:
	Gamma<6> gamma;
	static constexpr int M{N}; // power series terms
	static constexpr int64_t scale{1ll<<N};
	int64_t binom[N+1]{};
};
